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Abstract 

We develop a numerical algorithm for computing the effective drift and dif- 
fusivity of the steady-state behavior of an overdamped particle driven by a 
periodic potential whose amplitude is modulated in time by multiplicative 
noise and forced by additive Gaussian noise (the mathematical structure of 
a flashing Brownian motor). The numerical algorithm is based on a spectral 
decomposition of the solution to the Fokker-Planck equation with periodic 
boundary conditions and the cell problem which result from homogenization 
theory. We also show that the numerical method of Wang, Peskin, Elston 
(WPE, 2003) for computing said quantities is equivalent to that resulting 
from homogenization theory. We show how to adapt the WPE numerical 
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method to this problem by means of discretizing the multiphcative noise via 
a finite- volume method into a discrete-state Markov jump process which pre- 
serves many important properties of the original continuous-state process, 
such as its invariant distribution and detailed balance. Our numerical exper- 
iments show the effectiveness of both methods, and that the spectral method 
can have some efficiency advantages when treating multiplicative random 
noise, particularly with strong volatility. 
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Hermite polynomials 
2000 MSC: 60J25 
2000 MSC: 70K60 



1. Introduction 

We will develop and discuss numerical approaches to computing the long- 
time effective dynamics of a particle undergoing overdamped dynamics in a 
periodic potential, randomly modulated as a function of time, and driven 
additionally by thermal fluctuations. After rescaling the spatial variable 
with respect to the period L of the potential and temporal variable with 
respect to the time which the particle takes to fall from a maximum near 
to the minimum of the potential under the zero temperature dynamics, the 
equation of motion for such a particle can be written in one dimension in the 
form 

dX(t) = -0'(X(t))F(t)dt + v^dW^(t), (1) 

where X{t) denotes the particle position as a function of time, 9 = ksT /cf) is 
the ratio of the thermal energy ksT to the amplitude of potential variations 
0, is the periodic potential structure rescaled to have an order unity 
scale of variation, W{t) is a standard Brownian motion with {(iW{t)) = 
and {dW{t)dW{t')) = 5{t - t')dtdt', and F{t) describes the random tem- 
poral modulations of the potential, which may be assumed without loss of 
generality to have its amplitude (in mean and/or variance) normalized as 
desired. 

Equation ([1]) is an example of a Brownian motor , a class of stochastic 
systems which are used to characterize and analyze the mechanisms behind 
the functioning of biological molecular motors as well as to design artificial 
microscale and nanoscale machines j3,[5|. An overdamped dynamical descrip- 
tion (without inertia) is appropriate because of the small length scales and 
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physical values of the other parameters. The periodic environment reflects 
the ordered assembly of an extended microtubule, actin flber, or artiflcial 
substrate. Thermal fluctuations play an important but not entirely domi- 
nant role (so that 9 in practice tends to be somewhat but not much less than 
1; again a property of the microscale), and temporal modulations in the po- 
tential are induced by some external means (which typically does work on 
the particle). In biological settings, the potential in question is the binding 
potential between the molecular motor and the microtubule or actin flbers, 
and is modulated by chemical processes such as the binding of ATP or re- 
lease of phosphate, as well as physical processes such as an unbound "head" 
of the motor "searching" for and landing on a binding site (sj] . These physical 
and chemical processes proceed with effectively random delays because they 
typically rely on some component fluctuating under thermal effects until it 
manages to achieve a certain state so that the process goes forward. Con- 
sequently, continuous-time Markov chains (or sometimes renewal processes 
with more general waiting distributions 0]) are often used to describe the 
modulations F(t) j8-ll|, with each state corresponding to a possible geomet- 



ric conformation of the motor (except for its center of mass position, encoded 
by X{t)). In synthetic motors, the modulation F{t) may often be periodic 
by design 0, lij. We will here particularly focus on the special case of a 
continuous modulation by an Ornstein-Uhlenbeck process. Such continuous 
stochastic modulations of the potential could, for example, represent inter- 
ference of the motion of the molecular motor due to other organelles and 
structures in the cellular environment. 

A central practical question in the theory of Brownian motors is the 
overall long-time behavior of the particle. The periodicity of the potential 
and statistical homogeneity (or periodicity) of its modulations imply, through 



a central limit theorem argument 12|], that the statistics of X{t) at long time 



are Gaussian and characterized completely by the mean drift 



and diffusivity 



U.Umm (2) 

t-^oo t 

where (■) denotes a statistical average over all randomness. Of particular in- 
terest is how these transport parameters, characterizing the large-scale, long- 
time behavior of the motor, are related to the microscopic design parameters 
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(such as 9 and the structure of the potential and parameters characterizing 
the fluctuation F{t)) in the detailed stochastic differential equation model 
(II]). Analytical approaches are generally only possible in asymptotic limits 
(such as adiabatically slow or rapidly fluctuating modulations F{t)) ji], 
Some work has pursued such questions through direct Monte Carlo simula- 



tions [10|, of the stochastic differential equation ([T]), but this approach 
is rather expensive because the trajectories must be followed through many 
spatial periods (and typically also several realizations) and the nature of the 
Brownian motor (in particular 6' < 1) is such that the particle takes a sub- 



stantial amount of time to hop from one spatial period to another [15|, |16 
Accurate computations are further hampered by the slow convergence of a 
Monte Carlo simulation with respect to computational effort (square root ac- 
curacy gains with respect to simulation time and/or number of realizations). 

Deterministic numerical approaches can be alternatively developed based 
on the equivalence of the stochastic differential equation ^ for trajectories 
and the Fokker-Planck partial differential equation 

dtp{x, f, t) = -4 {-4>'{x)fp{x, f, t)) + 9d,,p{x, f, t) + C}p{x, f, t), (4) 

for the probability density p{x,f,t) of the particle position x at time t. In 
Eq. Cf is the infinitesmal generator operator associated to the Markov 
process F{t), and C*j its adjoint. Kostur ll5| developed finite element simula- 



tions of this equation with adaptive time stepping to achieve and estimate the 
long-time behavior for several canonical Brownian motor models, but notes 
that the periodicity creates some challenges for the implementation due to 
its spoiling of the band structure of the matrix formed by projecting the 
evolution operator in Eq. onto the finite element basis. Another tactic is 
to derive and numerically solve deterministic equations for the effective drift 
(|2]) and diffusivity ([3]) of the Brownian motor. Wang et al. [17'], Wang and 



Elston [18| designed an effective approach, which we will refer to as the WPE 
(Wang-Peskin-Elston) method and summarize in Section ^ for the case in 
which the potential modulations F{t) are governed by finite-state Markov 
chain dynamics. We in particular found this algorithm to be very efficient in 
mapping out the dependence of a two-state flashing ratchet model with re- 
spect to various underlying parameters [l[ . Another means of deriving direct 
deterministic equations for the effective drift and diffusivity is through ho- 



mogenization theory [19|, |20| . The resulting equations will be summarized in 



Section |3l Some relative virtues of this approach is that it can be developed 
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in a continuum framework, without committing to any particular discretiza- 
tion in advance, and follows a classical multiscale analysis. The equations 
in 17|,|l8|, on the other hand, are obtained after a particular numerically suit- 
able spatial discretization which together with the finite-state Markov chain 
structure of the modulations, induce a grand Markov chain structure to the 
dynamics. The derivation of the drift and effective diffusivity are obtained 
then by manipulations of the associated Kolmogorov equations featuring the 
transition rate matrices obtained by these discretizations. 

The formulas for the drift and the diffusion coefficient obtained using ho- 



mogenization theory can be rigorously justified [2l|, |22|. A natural question 



is whether other approaches that have been developed for the calculation 
of the drift and diffusion coefficients lead to formulas that are equivalent, 
at least in some appropriate asymptotic limit, to the ones obtained from 
homogenization theory. As examples we mention the calculation of the dif- 
fusion coefficient for a Brownian particle in a tilted per iodic potential using 
the mean first passage time (MFPT) approach |23l-l25l| and the calculation 
of transport coefficients (not only the diffusion coefficient) using the Green- 



Kubo theory [26|. The equivalence between the MFPT approach and the 



Green-Kubo theory with homogenization theory were investigated in [19 
and [27|, respectively. One of the goals of the present paper is to investigate 
the equivalence between the drift and diffusion coefficient formulas for the 
WPE algorithm and the (discretized) formulas derived from homogenization 
theory. The homogenization formulas can in particular be discretized in the 
same manner as WPE do at the beginning, with the result that the same 
discretized equation for the drift is obtained but different equations result for 
the effective diffusivity. After several numerical experiments verified that the 
two approaches achieved the same answer, we found that the WPE equations 
could in fact be derived by a variation of the homogenization argument by 
simply passing at one point to working with the adjoint of an equation. A 
unified framework capable of developing both the WPE and homogenization 
equations will be presented in Section HI 

Beyond simply providing another, possibly more transparent, framework 
for deriving the WPE equations, the homogenization approach affords some 
flexibility in the numerical discretization. The spatial discretization pur- 
sued by 171, |l8| is carefully designed to maintain the important property of 
detailed balance, and we do not seek to improve on this aspect. Rather, 
we consider how the effective drift and diffusivity for the Brownian motor 
equation ([T]) can be effectively computed when the temporal modulations of 
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the potential F{t) are Markovian and continuous in time. The prototypical 
example we shall examine is the Ornstein-Uhlenbeck process, which can be 
described equivalently as a Gaussian stationary random process with mean 
zero and correlation function 

{F{t')F{t' + t)) = ale-'/^ (5a) 

or as the solution of the stochastic differential equation 



dF{t) = --F{t)dt + J '^dWpit) (5b) 
r V T 

with -F(O) chosen as a mean zero Gaussian random variable with variance 
a p. In the above equations, Wpit) is another standard Brownian motion 
independent of W{t) in 1^, (Xp = (F(t)^) is the variance of F{t), and r is 
the correlation time of F{t). 

We will explore a spectral discretization of the state variable F which is 
closely related to the continued fraction method applied to many stochastic 



systems in Risken [28[ and to a neural network model in Acebron et al. [29 
but not, to our knowledge, to flashing ratchet equations ([1]). We compare 
in Section [5] the relative efficiency of the WPE and spectrally discretized ho- 
mogenization equations in computing the effective drift and diffusivity for the 
flashing ratchet ([T]) with modulations governed by the Ornstein-Uhlenbeck 
process (I5b[) . using Monte Carlo simulations as a point of reference for accu- 
racy. This model has been previously studied in the literature for the rapid 
decorrelation limit r J, and adiabatic limit r — )■ oo (see ^ and references 
therein) but, to our knowledge, the systematic computation of the effective 
diffusivity D is new. We will in particular study how the transport properties 
of a flashing ratchet with continuous Ornstein-Uhlenbeck modulations fl5bp 
compares with that of two-state flashing ratchets with the same correlation 
time and variance. 

We note finally that the methods developed here for the one-dimensional 
Brownian ratchet equation ([T]) can be generalized in principle to multiple 
dimensions, with a somewhat greater computational expense and more com- 
plex indexing of tensor products of Hermite polynomials. 

2. The Wang-Peskin-Elston Numerical Algorithm 

We begin by describing how the ideas from the Wang-Peskin-Elston (WPE) 



method [17|, ll8| can be adapted in order to compute the effective drift and 
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diffusivity of the flashing Brownian ratchet ([T]) modulated by the contin- 
uous Markov process F(t). First, F(t) is approximated by a finite state, 
continuous-time Markov chain F^{t) with state space {/nl^^^tt and transi- 
tion rate matrix K satisfying the property that all the row sums are zero 
^^,ggtt Knn' = and all non-diagonal entries are nonnegative. In the orig- 
inal formulation of the WPE method, F{t) is assumed to already be such 
a discrete-state Markov chain. The additive inverse of the negative diago- 
nal entries, —kn, defines the rate of leaving state n (inverse of the expected 
occupancy time), and the non- negative off-diagonal entries Knn' define the 
prochvity of jumping from state n to state n' in the sense that Knn'/kn de- 
fines the probability of such a jump whenever the Markov chain leaves state 



n 



30[. This step of discretizing the continuous Markov process F(t) must 
be performed carefully and in Subsection 15.11 we will elaborate on how we 
do this. Of course, if F{t) is already a finite-state Markov chain, this step is 
trivial. Next, the spatial variable x is discretized so that the Markov process 
{X{t),F{t)) defined in Section [1] is approximated by an extended Markov 
chain F'^{t)) on a finite state space (the Cartesian product of the dis- 

cretized state space of X and F). A key element to the WPE framework, 
particularly when employed for the purpose of trajectory simulation (not our 
focus here) in non-smooth (i.e., sawtooth) potentials, is the definition of the 
discretized Markov chain dynamics so that it preserves the detailed balance 
properties of the original equation ([1]). To explain this, we begin with the 
stochastic differential equation ([T]) with the Markov process F(t) replaced by 
a suitable Markov chain approximation F'^{t): 

dX{t) = -(f)'{X{t))F\t)dt + VWdWit). (6) 

We define p^{x,t) as the probability density for the semidiscretized process 
(X(t), F\t)),neSl, so that for any Borel set 5 G M and n G 

Prob {X{t) G 5, F\t) =n}= [ p"(x, t) dx. 

Jb 

The Fokker-Planck equation describing its evolution is given by 

^^^^^ = {<P\x)rP^{x, t) + ed,p^{x, t)) - Kpn{x, t)+J2 Kn'nP'^'ix, t). 

Next X{t) is also approximated by a discrete-state, continuous time Markov 
chain X<'(t) on a regular spatial grid Xij = j + iAx, i = 1, . . . , M^,, j G Z, 
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where we are assuming a spatial period of 1 and Aa; = 1 / ■ In this represen- 
tation of the grid, the parameter j indexes the real line by cells of length one 
(the normalized period of the potential 0,) while i indexes the grid points 
with separation distance l/M^ within each cell. The spatially discretized 
version of p"'{x,t) is indexed in a somewhat unorthodox way, convenient for 
the following developments. t) is defined to be the probability that the 
discretized joint process {X^{t), F'^(t)) takes values {xij, f^) at time t, and 
may also be interpreted as the probability that the semidiscretized process 
{X{t), F^{t)) takes values in the set (x^ — |Ax,Xjj + |Aa;] x {n}, and is 
therefore related to the probability density p"'{x,t) as: 

P7ij,t)= / p^{x,t)dx (7) 

This probability distribution is then represented in vectorial form p"(j, t) — 
{PiU: 0)P2 (i) t): - ■ ■ tPm^U: 0)- The discretized equations then read 

L-p^ij, t) + L:;p"(j -l,t) + L"_p-{j + 1, t) (8) 
+ ^ir„,„p"'(i,i), 

where the matrices appearing in this equation are given by 

[L% ^ - ^Fl^_^^^^ + Bly, + kn) forz = l,...,M., 

[Ll-M fori = l,...,M., 

[Ll+M =^+1/2' fori = l,...,M., 

(9) 

[1%, =0 for|i-z'|>2, 
['-+]i,M. =^M.+i/2, zero else, 

[^-]m^,i =^V2, zero else. 

Intuitively, the matrix L" is a discretization of the Fokker-Planck operator, 
while the matrices L" and L", acting upon the probability vector p^, are a 
discretization of the probability flux at the boundaries. The terms F^,^,^ and 
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-^1+1/2 represent the transition rates between adjacent cells, and are chosen 
such that for any frozen choice of F", the jump process X^, with master 
equation as given by Eq. ([8]) without the Kn'n term and the fc„ terms in 
L", satisfies the detailed balance condition with respect to the Boltzmann 
distribution e'^^^"^^^ . We refer the reader to Wang et al. and Wang and 



Elston [18| for further details of the numerical method. 

From the vectors and matrices defined above, we define supervectors and 
supermatrices with indices 1, . . . , M.w = x Np, where Np = \Sp\, so that 
the equation ([8]) can be expressed in the following abstract form: 

= Lp(j,t) + LMj - 1,^) + L-P(J + 1,^), 

More precisely, p"'{j,t) = [p]nM^+j i^)^ where [v]^ is used to denote the ith 
component of a supervector v for later convenience. Other supermatrices 
and supervectors are indexed similarly, with L a supermatrix representing 
the first and fourth terms in Eq. ([H]) (that is, the dynamics acting within a 

spatial period), while L+ and I are supermatrices representing, respectively, 

the second and third terms in Eq. (|8]) . These last two supermatrices will have 
a block diagonal form since they do not couple across different modulation 
states. 



Wang et al. [17|, Wang and Elston [18| show that the effective drift and 
diffusivity are obtained as the unique solutions to the following equations: 

Mw 

U = El(L+-L-)P^L' (10) 

1=1 

Mp' = 0, M = L + L_ + L+, (11) 



M 



w 



satisfying the normalization condition [p^]^ = 1, where A^iy is the total 

i=l 

number of discrete states. For the effective diffusivity, one must solve, 

D = -5^[(L+ + L_)p^ + 2(L+-L_)rL, (12a) 

i=l 

Mr = Up" - (L+ - L_)p^ (12b) 
with the normalization condition 

Mw 



i=l 
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3. Homogenization equations for effective drift and diffusion 

An application of the homogenization formalism from Pavliotis loj yields 



the following system of equations for the effective drift and diffusivity in the 
continuously modulated flashing ratchet model ([T]). We will not repeat the 
original derivation here since a unified derivation of the WPE and homoge- 
nization equations will be provided in Section |H 

First we define the infinitesimal generator of the Markov process {X{t), F{t)), 
to be understood as operating on functions defined on Sx x Sp where Sx is 
the spatial domain, here taken as the unit interval with periodic boundary 
conditions, and Sp is the state space of the potential modulations F{t). Pa- 
rameterizing the domain Sx by the variable x and Sp by the variable /, the 
infinitesimal generator reads: 

c = -<j)'{x)fd, + ed,, + Cf 

where Cf is the infinitesimal generator associated with the Markov process 
F{t). For the Ornstein-Uhlenbeck process fl5bl) . the associated state space is 
Sp = M. and the associated infinitesimal generator is: 

^f = l i-fdf + aldfj) . (14) 

To compute the effective drift and diffusivity, we first solve the stationary 
Fokker-Planck equation 

C*p{x,f) = 0, p{x + l,f) = p{x,f), (15) 

where 

C* = {<P\x)f-) + ed^^ + - {df if-) + aldff) 

T 

is the Fokker-Planck operator, defined as the adjoint of the infinitesimal 
generator. The stationary solution p{x, f) is to have periodic boundary con- 
ditions in X G [0, 1] and to satisfy the normalization condition: 

/»1 /"OO 

/ / p(x,/)d/da; = l. (16) 

Jo J -oo 

p{x, f) is here actually a reduced probability density 0, Sec. 2.4] associated 
with the Markov process {X{t), F{t)) in that only the relative position of 
X{t) with respect to the periodic potential is described; the information 
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concerning the nearest integer to X{t) is contracted out. This contraction is 
necessary for a nontrivial stationary probabihty distribution to be defined. 

Once p is computed, the effective drift is obtained rather simply as the 
drift coefficient in ([1]) averaged over spatial position with respect to the sta- 
tionary probability distribution: 

U = {-<P'{x)f), (17) 

where: 

p1 poo 

{9)p^ / / gixj)p{xj)dfdx. (18) 

Jo J-oo 

Next the following cell problem, which has the form of a Poisson equation, 
must be solved 

- Cx{x, f ) = -4>\x)f - U, x{x + 1, /) = x{x, /). (19) 

X must satisfy also the condition that it grows sufficiently slowly with respect 
to / so that: 

{\x{x,f)\^)p < oo. 

Moreover, the operator C has a one dimensional kernel of constants; the 
equation ( fT9l) does satisfy the solvability condition, and we impose an extra 
condition to fix a unique solution: 

= 0. (20) 

From here, the effective diffusivity D is computed as j3l| 

D = + {i-<P'ix)f - U) x)p + 2e{d.,x)p- (21) 

Note that the choice fl20l) does not affect the value of this formula for D, as 
it is invariant under the addition of constants to x- 

3.1. Comparison with Wang- Peskin-Elston framework 

We observe first of all some direct similarities between the homogeniza- 
tion equations just derived and those characterizing the Wang-Peskin-Elston 
framework summarized in Section [2j First of all, the equation Mp^ = 0, to- 
gether with its normalization can be readily understood as a finite volume 
discretization of the stationary Fokker-Planck equation f|T5|) with normaliza- 
tion f|T6l) . and describes exactly the reduced stationary probability distri- 
bution of the discretized system (X'*(t), F^{t)). The formula for the effective 
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drift ( ITU]) is also a simple discretization of equation ( 1T7|) . On the other hand, 
the formula ( I12ap for the effective diffusivity in the WPE approach does not 
appear to be a discretized approximation of the equation f l2T|) for the effec- 
tive diffusivity in the homogenization framework. In particular, the matrix 
M in ffTTl) and fll2bp is a discretization of the Fokker-Planck operator, while 



in order to find D from fl?I]) one must solve a problem involving the adjoint 
of this operator, as in f lTU]) . In Section HJ however, we will present a unified 
derivation of both the homogenization and WPE equations which explains 
the consistency of the results obtained with either method (up to numerical 
errors incurred by choice of discretization). 

For now, we proceed by noting that the homogenization equations have 
been derived without any prior discretization, so we may solve the equations 
for the effective drift and diffusivity using any method we please. We will 
next describe a spectral approach similar in spirit to the continued-fraction 



method of [28| developed in other contexts. 



3.2. Spectral decomposition 

We now develop a numerical algorithm to solve equations (fT5|) and (1191) 
through a spectral decomposition in terms of Hermite polynomials (in /) and 
Fourier series (in x) of their respective solutions p{x,f) and x{^if)- This 



method is analogous to that presented in [28| for computing the effective 
drift of a particle on a tilted periodic potential, including its extension to 
computing effective diffusivity ISl], and which was employed in Latorre et al. 
32| for numerical comparison against theoretical expansions of the trans- 
port coefficients with respect to the strength of the tilt and corresponding 
corrections to the Einstein relation. We commence by writing equation f[T5|) 
as: 

{(j)'{x)fp{x, /)) + ed^^pix, f) + -C}p{x, /) = 0, 

where £j is the adjoint operator of £/ given in Eq. (IT^ . The invariant 
probability density of the F{t) dynamics is given by the solution of 

qPFif) = 0, 

which is readily seen to be prif) = {27Tap)^^^'^e~^'^^^'^'^'f'\ We write the 
reduced stationary distribution for the dynamics {X{t), F{t)) by factoring 
out this expression 

p{xj) = pF{f)'rc{xJ), 
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where n{x, f) is just the conditional (reduced probabihty) density of X{t) 
given F(t) = f. After substituting this ansatz in equation (IT^ we are left 
with 

(0'(x)/7r(x, /) + eOM^, /)) + Um^, f) = 0. (22) 

A convenient basis which diagonalizes the £/ operator (for Ornstein-Uhlenbeck 
dynamics ( I5b[) ) and is orthonormal with respect to pi? is given by the Hermite 
polynomials {Hn{f)}'^=Q, which can be defined through 

and have the properties: 

CpHnif) = -nH„if), n = 0,l,2,... 

/oo 
H^if)HUf)PF{f)df = 5nm. 
-oo 

We now expand the solution n{x,y) to the equation ( l22i) with respect to 
these Hermite polynomials: 

oo 
n=0 

This representation of 7t{x, /) decomposes equation fl22]) into an infinite sys- 
tem of coupled ordinary differential equations, 

apd^r {(f)'{x)TTi{x)) + ed^^TToix) = 0, (23a) 

£^7r„+i(a;) + £„7r„(a;) + £+7r„_i(a;) = 0, n = 1, 2, . . . , (23b) 

with 

£„7r„+i(a;) = \/ {n + l)aFd^ (0'(x)7r„+i(x)) , 
C+Tin-iix) = ^/naF^^{(f)'{x)^Tn-l{x)) , (24) 
CnTCnix) = [dd^^ - TiT^^) 7r„(a;). 
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Since iin{x,f) must be periodic functions, an obvious choice for solving the 
above infinite system of equations is to express vr„(x) in terms of Fourier 
modes, 



j=-oo 



(25) 



This representation leads to an algebraic system which is solved by truncating 
the series. 



Ns Ms 



n=0 j=-Ms 



We refer the reader to Appendix A for the details of the algorithm and how 
the system is solved. In the end, U is computed from Equation f|T7|) as 



D is computed in a similar way |31i; we expand the solution xi^if) with 
respect to the Hermite polynomials in / and Fourier series in x, and truncate: 

Ns Ms 
n=0 j=-Ms 

After solving for these expansion coefficients through a projection of the 
governing equations (|T5l) and f|T9|) . we evaluate the effective diffusivity as 



D 



N 



n=0 

N M 



Ms 



Ej j+1 j j—1 I j j'+l j j—1 



.3=- Ms 



(26) 



71=0 j=-M 



Details of the derivation can be found in Appendix A 



4. Equivalence between the WPE Numerical Algorithm and Ho- 
mogenization Theory 

We show in this section how the WPE algorithm described in Section [2] 
and the homogenization theory in Section |3] can be obtained through a unified 
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multiscale derivation, in which one simply chooses at one point between 
working with an equation or its adjoint, and then of course a specific choice 
of discretization. 

We begin by considering the fiashing ratchet model ([1]) with arbitrary 
Markov process F{t) modulating the potential, and do not yet impose any 
discretization. The derivation of the homogenized expressions f lT7|) and f l21l) 



for the effective transport coefficients in Pavliotis [19| pursued a multiple 
scale analysis of the backward-Kolmogorov equation, 

du{x,f,t) 

— = Cu{xJ,t). 

In order to establish the equivalence between the WPE numerical algorithm 
and the homogenization theory approach, we find it more convenient to in- 
stead apply the multiscale technique to the forward-Kolmogorov or Fokker- 
Planck equation associated with ([1]), namely 

dtp{x, /, t) = C*p{x, /, t) = (0'(x)/p(-)) + Od^^p{x, f, t) + £}p(x, /, t), 

(27) 

where C*j is the adjoint of the infinitesimal generator of F with state space 

4-1- Equivalence of Drift Formulas 

We begin by seeking a coarse-grained description on advectively-rescaled 
large space and time scales 

PA{xJ,t) = e~'p{x/eJ,t/e) (28) 

and seek a two-space scale solution of the form 

PA(x,/,t) = pSs,A(^,e,/,t) , (29) 

with small parameter < e ^ 1 denoting the separation of scales between 
the coarse-grained observation scale and the period, ^ the small-scale space 
variable. We don't include a small-scale time variable because the structure 
of the dynamics is such that the statistical distribution should approach a 
quasi-steady state on the small-scales 0, Sec. 2.4], and we are not inter- 
ested in resolving the transient evolution of the small-scales from the initial 
data. We seek solutions Pmsa which have periodicity Pmsa(^'^ + l;/;^) = 
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Pms a(^'^' /, t) in the small-scale variable ^ corresponding to the periodicity 
of the potential. With the chain rule, we find that a solution of 



A 



( 0'(O/p£kA ) + (^d^^A + ^}P)^SA 



+ 



(30) 



yields, through Eqs. f l28|) and (|29|1 . a solution to Eq. ( l27|l . Substituting next 
a perturbation expansion 

Pms a(^' ^' /' ^) = PAo{x, ^, /, t) + epAi{x, ^J,t) + ... 



(31a) 
(31b) 

(31c) 



into Eq. (!30|) . we obtain the asymptotic hierarchy 

O(e-i) : = CIpao, 
0(1) : 9tPA0 = QpAi + (f)'{Of9xPAo + 2ed^^pAi, 



where we have defined the fundamental operator on small-scale variables: 



In solving these equations, we use the following solvability condition [33[: 
The equation 

>CS^(e,/) = Me,/), (32) 
has a periodic solution g{C,, f) = g{C,+^, f) only when the solvability condition 



r / Me,/)d/de=o 

Jo JSf 



is satisfied. When this condition holds, Eq. fl32|) has a one-parameter family 
of solutions g{^, f) = gp{^, f) + CTTo{^, /) where c is an arbitrary real constant, 
and vro(^, /) is defined as the unique real, periodic solution 7ro(C, /) = 7i^o(C + 
1, /) of the homogenous equation 



with the normalization 



= 



7ro(e,/)d/de = l. 



(33a) 



(33b) 
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In the case in which F has discrete state space 5*/, the integral over Sf 
should be replaced by a sum over states. This solvability condition is de- 
rived from the fact that the operator Cq is elliptic, with one-dimensional 
null space spanned by ttq and one-dimensional adjoint null space spanned by 
constants |33|. 

Applying this solvability condition to the 0(e^^) equation in Eq. f l?!]) 
yields the result that Pao{.x,^, f ,t) = ttq^C^, f)c{x,t) for some function of 
large-scale variables, c{x,t), to be determined. Substituting this expression 
for pao into the 0(1) equation in Eq. ( 13T|) . and then imposing the solvability 
condition produces the result that 

%^ = V.(Ue(.,t)) 
which is of course a simple advection equation with drift velocity 

U = r / -0'(O/vro(e,/)d/de 

Jo J Sf 

This recovers the homogenization formula (ITTj) for the effective drift of 
the motor. We show now how upon an appropriate discretization of physical 
and state space, the WPE formulas (fTTj) and (fTOl) can be recovered. First 
of all, the equation ( ITTl) is formally just a discretization of the equation (1151) 
for the stationary distribution of the position of the motor. To make more 
precise contact with the choice of supermatrices L, L_|_, and I used in Wang 



et al. [17|, Wang and Elston [18[, we observe that equation (I33ap can be 
written as. 



+ £}7ro(^,/) = 0, (34) 



A direct finite- volume numerical discretization leads to a numerical scheme 
of the type developed in js^l , and a further consistent approximation of the 
coefficients in the resulting algebraic equations leads to the equation ( ITTll 
with the forward and backward rates in (|9]) equal to those in Wang et al. 



17j . Wang and Elston 18|]. To derive the WPE formula f lTOj) for the effective 
drift, we integrate Eq. over the noise variable /, noting the integral over 
the £j term vanishes due to integration of a derivative of a function with 
decaying values as |/| — )■ oo, to obtain the following relation: 

d 



^e-MO/^9,(e/^«)/^o(e,/))d/ = 0, (35) 

Sf 
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This implies that the integral, which is nothing but the averaged net spatial 
flux induced by the stationary distribution at a position ^, is constant on 
the cell [0,1], a statement which can alternatively be derived by physical 
considerations We next observe that he expression f|T7|) for U can be 
written as 

Jo JSf 

which, upon integration by parts and the observation that all factors in the 
integrand are periodic, leads to 

U = -/' / ^e"/^«)/^9,(e/^«)/%(e,/))d/de 

^0 JSf 

But by Eq. f l35|) . the integrand of the ^ integral in this expression is indepen- 
dent of ^, and therefore everywhere equal to its value at the period boundary 
^ = 1, so we may equivalently write 



U = - 



Sf 



(36) 



In continuous variables, this is precisely the flux of probability through the 
period boundary at ^ = 1, which again has a clear physical interpreta- 
tion The equation (fTOj) can be similarly understood, through the general 

expression of the supermatrices L+ and I in terms of forward and backward 

transition rates, as the discretized form of the net probability flux through 
the period boundary at ^ = 1. More precisely, we approximate the deriva- 
tive in Eq. [36] through a standard centered flnite-difference at the points 
^ = 1 -|- Ax/ 2 and ^ = 1 — Ax/2, and discretize the state space Sp as 
described at the beginning of Section |2] to obtain an expression of the form 

E KW)m.-^Up1i] (37) 

for suitable coefficients A'^^^ and We have used the periodicity of the 

stationary distribution to identify 7ro(l + Ax/2,-) = 7ro(Ax/2,-), which is 
then approximated in standard fashion (Eq. ([7])) in terms of (temporarily 
relaxing supervector indexing for clarity). As before, a further consistent 
approximation of the coefficients A^^^ and Al^_ leads to the WPE formula f lTO]) 
for U, with the supermatrix coefficients defined precisely as in Wang et al. 



17f . Wang and Elston [18 
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4-2. Equivalence of Diffusivity Formulas 

We turn next to the derivation of the effective diffusivity D for homog- 
enization theory ( !T2|) and the WPE equivalent (fT2|) . which unhke the drift, 
appears not to be a simple discretization of the homogenization formula. To 
this end, we rescale diffusively to large time and space scales, centered about 
the net drifting motion: 

We seek a solution of the form 

pD{xJ,t) = p'^^^ix^^J^t) 

5=(x+Ui)/e 

periodic with respect to the small space variable ^. Note we have defined 
the small scale variable ^ to undo the Galilean transformation to avoid the 
need to add a fast time scale arising from trivial advection of the stationary 
small-scale structure. By the chain rule, we can generate suitable solutions 
of Eq. (127|) through periodic solutions to the multiscale transformed Fokker- 
Planck equation: 



5«((0'(O/p£kD)+^%PtkD 

^-1 



(38) 



Upon substituting the perturbation expansion 

PmsA^^ t) = PDo{x, ^, /, t) + epDi{x, ^, /, t) + e^pD2{x, ^,f,t) + ... (39) 
into Eq. fl55]) . we find by equating equal powers of e, 
0(e-2) : = C*,pDo, 



0(6-1) 
0(1) 



= C*QpDi + CIpdi, 



(40) 
(41) 
(42) 



where. 



n 



:9a0'(O/-) + ^% + ^/, 
jCI =9,((0'(O/ + u) ■) + 2ed^x. 

Using the same solvability condition as above, the 0(e~^) equation implies 



that 



Pdo{x, ^, /, t) = c(x, t)7ro(^, /), 



(43) 
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where c(-) is a function to be determined, and ttq is defined as in Eq. 
Equation (HT]) reads, 

- QpDi = C\Pm = WiOf + U) TTo + 2ed^T^o] dxc, (44) 
which is automatically solvable since the right hand side satisfies 



/ CIpq df d^ = d^c 

JSf \JSf 

= 0, 



1 








5=o_ 





because of the periodicity of vro(^, /) and the definition ( 1T7|) of U. As the 
variables x and t enter as parameters in (jH]), we can express the solution in 
the form 

PDiix, /, t) = $(e, f)d^cix, t) + 7ro(e, /)ci(x, t), (45) 

where Ci{x,t) is for now an arbitrary function and $(z, /) is the unique 
solution of the equation 

- Cm^, f) = (0'(O/ + U) 7ro(e, /) + 290^^, f) (46) 
which is periodic, satisfies (|$p)p < oo, and with integral chosen to be: 

f [ $(^,/)d/de = - r / e7ro(e,/)d/de. (47) 

Jo JSf Jo JSf 

We note that any constant could have been chosen on the right hand side, 
without affecting the subsequent derivation of the formula (HSj) for the effec- 
tive diffusivity, but our particular choice will facilitate connection with the 
WPE formula (fT2|). 

This is an adjoint equivalent of the cell problem f lT9|) that arises from 
the homogenization analysis of the backward-Kolmogorov equation. Upon 
substituting the results fH5]) and f HS]) into Eq. f H2]) . we obtain 

7ro(e, mc{x, t) = {9noi^, f) + (0'(O/ + U) $(e, /) + 290^^^, /))Cc(x, t) 

+ ( + U) 7ro(e, /) + 29d^7ro{^, /))9,Ci(x, t) + CIPd2{-)- 

The solvability condition for this equation then implies that c{x, t) satisfies 
the diffusion equation 

dtc{x,t) = Ddl^c{x,t), 
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where D is given by 



D= r/ (^vro(e,/) + (0'(O/ + U)$(e,/) + 2ea^$(e,/))d/de 
Jo JSf 

= e+ f [ (</>'(o/ + u)$(e,/)d/de, (48) 

Jo JSf 



where we have used the normahzation fl33bp of ttq and the periodicity of 



$ in the last equahty. This is a somewhat different expression for D as 
was obtained in Eq. (!2T|) from the same multiscale technique apphed to the 



backward Kolmogorov equation [19[. A rigorous justification of the above 
formal multiple scales analysis, together with a proof of the fact that the 
diffusion coefficient is finite, can be obtained using tools from stochastic 
analysis, such as the martingale central limit theorem, together with a careful 



study of the cell problem [35 



Before showing how the expression f H8|) for the diffusivity is related to 
the WPE algorithm, we prove directly the equivalence with the equation ( 12T|) 



that arise from the original homogenization theory from Pavliotis [19[. The 
latter are expressed in terms of an auxiliary field xiC^ /)? which satisfies the 
cell problem 

with periodicity in ^ and integrability (|xP)p < oo. 
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From equation (HSi) we have, 

D =0 + [' [ (0'(O/ + U)$(^,/)d/de, 

^0 JSf 

= 9+ [' [ $(e,/)Ax(e,/)d/de, 

^0 JSf 



= 9+ [' [ x{^,mm^,f)dfd^, 

Jo JSf 

from iS]) =e + /) + U)7ro(e, /) + 268^^^ /)) d/ d^, 

= ^ + I' I - U)x(e, /)vro(e, /)) - 2^^x(e, /)5cvro(e, /) d/ d^, 

= ^ + l' I {(-^'(Of - U)x(e, /)7ro(e, /)) + 2^vro(e, /) d/ d^, 

= 9 + - U)x).o + 2^(%).o- 

The above expression is precisely equation ( 12T|) for D found via the original 



homogenization theory in [19 



The equivalence of the expression ( l48ll for the effective diffusivity with 
that of the WPE numerical method can be established as follows. From 
equation ( H6l) . 



-Cm^, f) = + U) 7ro(e, /) + 29d^7ro{^, /). 

Now define 

i?(e,/) = -($(e,/) + evro(e,/)). (49) 

It is easy to verify that R satisfies the boundary condition 

i?(e + l,/) = i?(e,/)-vro(e,/) (50) 

and the equation, 

£*i?(e,/) = U7ro, (51) 
and moreover from Eqs. ( H7|) and ffT6|) . that 



r / i?(e, /)d/de = 0. 

^0 JSf 



(52) 
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The same argument based on finite-volume discretizations used to connect 
the continuous equation fl33ap with the discretized WPE equation (ITU]) , show 
that fll2bp is simply a discretization of Eq. (|5T1) . The matrix M is a discretiza- 
tion of the operator £q corresponding to periodic boundary conditions, and 
r = i?2, • • • , Rmw)^ ^ P' = (^1' 7r2, . . . , T^Mwf are the corresponding dis- 
cretized approximations of R{^, f) and vro(^, /). The term — (L_|_ — I )p^ is 

a correction term arising from the fact that R{^, /) does not satisfy periodic 
boundary conditions, but rather Eq. (l50l) . Thus, the discretized terms corre- 
sponding to "fluxes of R" from outside the period domain [0, 1] must involve 
this shift of ivTo relative to the case in which R is periodic and fluxes in from 
the left/right are equated to fluxes out of the right/left boundary. Also, the 
normalization condition (1131) is clearly a direct discretization of the integral 
condition fl52|) . 

Finally, we show how f H8|) leads to the WPE formula fll2al) for the effective 
diffusivity. First, we note that by periodicity of $(^, /) and the integration 
over a complete spatial period in $, as well as Eq. fH9l) and the normaliza- 
tions ffTBl) and f[T^ . we can rewrite Eq. fHS]) as: 

B = 9+ f [ (0'(e)/ + ^95 + u)$(e,/)d/de 

Jo JSf 

= ^ - / V (0'(O/ + + U) m, f) + e^o(e, /)) d/de 

^0 JSf 

Jo JSf jo JSf 

- f I eno{^,f)dfd^-u\ [' [ R{^,f)dM+ f [ e^o(e,/)d/de 
Jo Jsf Uo Jsf Jo Jsf 



'Sf 

9[e-^^/'d,{ef^/'R{^,f))]dfd^ 



JSf 

\ I eyf^''d,{ef^'\,{ij))]dfdi 

Jsf 

f I ^vro(e,/)d/de (53) 
Jo JSf 

Next we integrate by parts in the first integral, noting the periodicity of all 
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factors in the integrand except i?, which satisfies Eq. ( 15U]) . to obtain: 



'Sf 



(e^^/^i?)] d/de 



Sf 



+ 



from (1511) 



+ 



Sf 
1 



Sf 

f 

Sf 



^d^ {9e- 
-f^"d^ (e 

eUTTod/de 



(54) 

'ae^^/^i?))d/de 
d/l 

5=1 



Noting from our argument from Eq. (|35l) . the integral over in the second 
term in the last expression of Eq. ( |53l) is in fact constant with respect to 
we can trivially integrate over ^ to obtain: 



-I q 

'0 JSf 



Sf 



^e-^<^/^95 (e^^/^o) df 
(55)" 

Combining then Eqs. ( 153|) . ( IMI) . and ( 155|) . we obtain the following expression 
for the effective diffusivity in terms of continuum variables which, analogously 
to the formula Eq. ( !36|) for the effective drift, only involves evaluations of 
"fluxes" at the period boundary rather than integration with respect to the 
spatial variable ^: 



D = - 



df 



€=1 



(56) 



Were (i? + Ittq) a periodic function, then by the same argument as above 
which interpreted the operator in Eq. ( 136|) as a net spatial flux, and the 
matrices L+ — I as a corresponding discretization, we would say that the 



right hand side of Eq. (1561) is just J2 



Mw 
1 



-L-)(r + ip^ 



But R is 



not periodic, and this argument is flawed because while L+r does serve as an 
appropriate discretization of the "rightward spatial flux" of R across ^ = 1, 

I r describes the "leftward spatial flux" of R across ^ = 0, which is not the 

same as the "leftward flux" of R across ^ = 1 due to the lack of periodicity 
of R. Rather, since by Eq. ([50D, R(l + ^J) = /) - vro(^, /), we should 
discretize the "leftward spatial flux" of R across = 1 as I (r — p^). Then 



€=1 
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taking the net spatial flux at ^ = 1 as the "rightward spatial flux" minus 
the "leftward spatial flux,"' integrated over the noise variable /, we would 
discretize Eq. fl56|) as: 



D 



Mw 



Ur-L_(r-p^) + (LH 



^ Mw 



L-)r+(L, 



L-)-p' 



in agreement with the WPE expression fll2al) . 

We can make this argument somewhat more concrete, as we did for U 
around the discussion of Eq. fl37|) . by approximating the last expression in 
Eq. f l56|) by a centered finite-difference using the points ^ = l + Ax/2 and ^ = 
1 - Ax/2, but now we must notice that R{1 + Ax/2) = R{Ax/2) - tiq{Ax/2) 
(while 7ro(l + Ax/2) = 7ro(Ax/2)), to obtain: 



D 



Mx/ n _|_ 1 (r,^\n 



nGS« 



where we temporarily suspend supervector indexing of r and p®. The same 
consistent approximation of the coefficients A^^J^: and Al^_ as in our discussion 
of the effective drift gives precisely the WPE formula fll2al) for the effective 
diffusivity, with supermatrix coefficients defined precisely as in Wang et al. 



17[, Wang and Elston [18| 



5. Numerical Results 

In this section we will explore the efficacy of the WPE and homogeniza- 
tion approaches in simulating a flashing ratchet ([6]) where the potential is 
modulated by an Ornstein-Uhlenbeck process (I5b[) . Though we have shown 
that the formulas for the effective drift and diffusivity are formally equiva- 
lent for the two methods, their implementations differ in their discretization. 
In particular, the homogenization algorithm will discretize both space and 
the random potential modulation through a spectral expansion, as discussed 
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in Section |21 On the other hand, the WPE algorithm will discretize both 
space and the random potential modulation through regularly spaced grids, 
in such a way that the stochastic differential system ([6]) is approximated by a 
finite-state Markov chain which preserves detailed balance (Subsection I5.ip . 
We also examine the theoretical question of how the transport properties of 
the motor compare under continuous-state or discrete-state potential modu- 
lations with equivalent low order statistics (Subsection 15. 2p . The numerical 
studies for both the comparison of the algorithms and the discrete-state and 
continuous-state modulations are presented together in Subsection 15. 3[ 

5.1. Discrete-state approximation of the Ornstein-Uhlenbeck process. 
The equation of motion of the flashing ratchet is given by, 

dX{t) = -(j)'{X{t))F{t)dt + V2edW{t), 

where F{t) is the external modulation, which we now fix as the Ornstein- 
Uhlenbeck (OU) process 



dF(t) = --F(t)dt + J^dWFit). 
T V r 

The means for computing the effective drift and diffusivity for this system us- 
ing homogenization theory, including their discrete, computable approxima- 
tion, was presented in Section [3l The WPE numerical method for computing 
the transport coefficients, on the other hand, is formulated in Wang et al. 



17l |. Wang and Elston [18] for flashing ratchets where F{t) is a continuous- 
time Markov chain with a finite state space. Consequently, to implement 
this approach on the continuously modulated model (jH]), we must somehow 
approximate the continuous dynamics of the OU-process by a finite-state, 
continuous-time Markov chain. This is usually done by approximating ei- 
ther the backward-Kolmogorov or forward-Kolmogorov equation with some 
finite-difference method. This approach, however, may give rise to some in- 
consistencies in the approximation. For instance, if the grid size is not small 
enough the jump rates may not be positive, and some important properties 
of the original continuous process, such as those concerning its invariant dis- 



tribution, may be lost. Recently, in [3J], a numerical method was developed 



in order to consistently approximate a continuous-state stochastic process 
with a Markov jump process by using a finite-volume method. Here, we 



In fact, the WPE numerical method can also be derived using this approach. 
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will present this approach for the simple one- dimensional process F{t). We 
begin by writing the backward-Kolmogorov equation for F{t), 

du{f,t) _ 1 fj^^ ^ ^2 



dt t\ df ^dp 

Defining (3 = and = /^/2, the backward-Kolmogorov equation can 
be rewritten as 

MM _ 4em/)^ fe-^^(^)#.(/,t)') . (57) 



dt T 9f \ df 

Once written in this form, a finite-volume method can be used to discretize 
the equation. In one dimension and for a uniform grid, a simple finite- 
difference scheme can be used to obtain the same approximation. This is 



presented in detail in [Appendix B[ In the end, the approximation of the 



backward-Kolmogorov equation can be expressed as, 

^u(t) = Lu(t), 

where the entries of the matrix L are given as, 

[L]n,n' = 



-{Kn,n+1 + Kn,n-l) H U = Tl' , 

Kn,n+i if n' = n + l, 

Kn^n-i if n' = n — 1, 

otherwise. 



where the Kny are nonnegative constants with expressions given in Eq. (IB. ID . 
and u(t) = {ui{t) , U2(t) , . . . ,uiyp(t)) are the pointwise approximation of the 
solution u{f, t). The spatial discretization for the WPE method then follows 



the standard procedure described in Wang et al. jl71, Wang and Elston [18 
building upon this finite-state Markov chain approximation for F{t), which 
we denote F^{t). We will refer to the resulting numerical method, extending 
the WPE ideas with the finite-volume discretization procedure of Latorre 



et al. j3J] to handle the discretization of the stochastic process F{t), will be 



referred to as the 'WPE-based" method in the following discussion. 

5.2. Comparison between discrete- state and continuous- state flashing ratchet. 

We next turn to the question of how much of a difference a continuous- 
state modulation of the flashing ratchet fl5b|) has on the transport properties 
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relative to a discrete-state dynamics with the same low-order statistics. Of 
course we expect that with sufficiently many discrete states, the transport 
properties should be relatively insensitive to the discretization, so we set the 
comparison most starkly by comparing the Ornstein-Uhlenbeck modulation 
with a dichotomous Markov-chain modulation Foit) taking values {/i,/2} 
with transition rates between the two states given by ki2 and ^21- We choose 
these parameters to mimic the Ornstein-Uhlenbeck process as closely as pos- 
sible. First, because the OU-process is symmetric about the origin, we set 
fi = — /2 = / and ki2 = /C21 = k. This makes the dichotomous process 
have mean zero, as does the OU-process. We next ask that both the discrete 
and continuous process have the same correlation function, assuming both 
are initialized with respect to their stationary distributions. The OU-process 



has correlation function 36 



{F{t')F{t' + t)) = ale-'/^ 
whereas the dichotomous Markov chain has correlation function 

{FD{t')FD{t' + t)) = pe-''\ 

We set then, 

f = ap, k = — . 

These restrictions completely determine the Markov chain Fd. 

5.3. Comparison of WPE-hased and Homogenization Algorithms for Contin- 
uous Potential Modulations 

We explore the performance of the homogenization algorithm and the 
WPE-based method for a rather simple example in which the potential is 
sinusoidal (j){x) = —-^ cos2'kx (for which the effective drift U = 0). In 
Figures da] and [lb] we present the results of the computations for the effective 
diffusivity as a function of the variance ap of the modulations F{t) for two 
different values of the correlation time r. For r = 0.01 (Figure [Ta|) we observe 
an enhancement of diffusivity (i.e., D > while for r = 10 (Figure [Tb]l we 
observe a suppression of diffusivity (i.e., D < 6). In these figures, we have 
used Ms = 20 (41 Fourier coefficients) and Ng = 30 Hermite polynomials for 
the spectral method, while using M^ = 500 grid points in the x-direction and 
Np = 21 grid points in the /-direction for the WPE-based method (resulting 
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in a Markov jump process with 21 states; see Appendix B for how A/ is 
chosen.) 

The Monte Carlo simulations were performed by an Euler-Maruyama 
discretization of the SDE ([6]), with a time step At = 0.001 and an ensemble 
average over 1000 independent simulations after a large number of time steps, 
which varies depending on the parameters of the simulation (see the figure 
captions for the actual number). The diffusivity D for the flashing ratchet 
with dichotomous noise was also computed via the WPE-based method, using 
the same number of grid points. 
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Figure 1: Effective diffusivity D as a function of multiplicative noise variance 
ap {6 = 0.1), computed for the OU- flashing ratchet IQ with the spectral 
homogenization algorithm (solid line), the finite- volume adaptation of the 
WPE-based numerical algorithm {Np = 21, dotted line.) The dash-dot line 
indicates the effective diffusivity for a flashing ratchet with dichotomous noise 
with same mean and correlation function as the OU- flashing ratchet. Monte 
Carlo simulations after ([HD 3 X 10^ (([Tb])-([Idl)) 10^ time steps (solid line 
with one standard deviation error bars). 



We can observe from these figures how the actual number of states for 
the multiplicative noise plays a fundamental role as the fluctuations of F 
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become larger, especially for larger values of r. In Figures [Tc] and [Td] we 
present computations of the effective diffusivity D for intermediate values of 
r, where this phenomenon is also observed. The parameters for the algorithm 
in the numerical simulations are the same as before. 

We can observe from the figures that each of the methods are computing 
the effective diffusivity consistently for the parameter ranges explored, and 
that the behavior of the motor particle is sensitive to whether the flashing 
ratchet is discrete or continuous precisely when the correlation time is large 
and the amplitude of the potential modulations is not small (in our rescaled 
units). 

For another perspective on the results, we study next the behavior of 
the effective diffusivity as a function of the parameters r, D, and 9, where 
D = (tI^t. This latter parameter characterizes the strength of the noise some- 
what differently than the simple amplitude by taking into account also the 
correlation time. D can be thought of as a crude scaling estimate, from ki- 
netic theory principles, of the enhancement of the diffusivity of the motor 
particle due to the flashing ratchet, and should be accurate (up to constant 
prefactor) for the case of low Kubo number 37| in which the decorrelation 
in the motion of the motor particle is determined essentially by the temporal 
decorrelation of the amplitude modulation F{t) rather than spatial decorre- 
lation through motion across the potential landscape 0(x). 

In Figure [2] we present our findings when we fix D = 1. The param- 
eters for the algorithm are the same as before. We see first of all that 
the homogenization algorithm remains in good agreement with the Monte 
Carlo simulations throughout the range of correlation times presented, and 
that the effective diffusivity so computed agrees with the intuition described 
above that D ~ CD for small correlation time r (with some order unity 
constant C ). On the other hand, the WPE-based algorithm with Np = 21 
fails to follow the Monte Carlo simulations when the correlation time r of the 
Gaussian noise is very small (and consequently the noise amplitude ap is very 
large). In this scenario one must increase the number of states in the Markov 
chain approximation of the OU-process to obtain accurate results with the 
WPE-based algorithm. We note also the related observation that for small 
correlation times r and fixed D = CpT, the behavior of the motor particle 
becomes very sensitive to whether the potential modulations are continuous 
or discrete. This regime corresponds to a limit in which F{t) approaches 
white noise with correlations {F{t')F{t + t')) = T)6{t). Combining the ob- 
servations from Figures [T] and [H we see that the dichotomous Markov chain 
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Figure 2: D as a function of r, with D = apT = 1, 6' = 0.1, computed 
for the OU-flashing ratchet (|6]) with the spectral homogenization algorithm 
(solid line), the WPE-based numerical algorithm {Np = 21, dotted line), and 
Monte Carlo simulations after 2 x 10^ time steps (solid line with one standard 
deviation error bars). The dash-dot line indicates the effective diffusivity for 
a flashing ratchet with dichotomous noise with the same mean and correlation 
function as the OU-fiashing ratchet. 



approximation to the continuous Ornstein-Uhlenback process for the random 
potential modulations creates similar behavior for the motor particle, except 
when the amplitude ap of the fluctuations is larger than some critical value 
which decreases with the correlation time r. 

5.4- Cost Comparison between the Spectral and WPE-based numerical Meth- 
ods 

We present now a comparison of how the solution of the numerical meth- 
ods presented above converge with respect to the number of elements taken 
in the approximation. As we saw in the previous section, the number of 
states Np in the discrete approximation of the OU process plays an im- 
portant role in the accuracy of the WPE-based method, especially for large 
values of a p. This should come as no surprise, for the approximation is based 
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on a finite-volume approximation of tlie backward-Kolmogorov equation of 
the OU process. Then the factor 1/Np is proportional to the grid size A/. 
We present then convergence comparisons between the two methods for re- 
finements along both the spatial and noise variable directions. To represent 
the spatial discretization, we use the number of grid elements, Mx-, taken in 
the WPE approximation for the X{t) process and the number of Fourier el- 
ements, Ms, taken in the truncation in the spectral algorithm. For the noise 
variable discretization, we represent the computational effort with the num- 
ber of grid points, Np, taken in the finite- volume approximation of the OU 
process and the number of Hermite polynomials, Ns-, taken in the spectral 
algorithm. In Figure [3] we present how the error in the numerical solution 
of the spectral algorithm is reduced as we increase the number of Fourier 
elements in the truncation while keeping the number of Hermite polynomials 
constant. The error is computed as usual as 

error(M,) = |DAw-D(M,)|, 

where DM^ax is the solution using a large number of Fourier elements Mmax (in 
this case it is double the number of the last simulation point), and D{Ms) is 
the solution computed using Mg Fourier terms (analogously, grid points 
in the x coordinate for the WPE-based method). In Figure S] we present 
the same experiment for the WPE-based method. In this case, we compute 
the error in the solution as we increase the number of grid points in 
the X direction, while keeping fixed the number of grid points Np in the 
F-direction. 

In Figures we perform a similar experiment but now increasing the 
number of Hermite polynomials Ng in the solution of the spectral method 
while keeping the number of Fourier elements fixed. Analogously, in Figure 
[6] we increase the number of grid elements Np (equivalently to decreasing the 
grid size A/) while keeping the grid size Ax constant. 

We can clearly see the Ax^-convergence in the WPE-based method (as 
well as A/^-convergence), which is characteristic of 2nd-order finite- volume 
approximations. On the other hand, it is clear how the spectral method con- 
verges faster as the number of spectral elements are increased. 

A natural question now is how the error in both methods converges as 
the cost of the numerical method is increased. Although a careful analysis of 
the numerical cost (as given by the number of flops, for instance) is beyond 
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Figure 3: Error in the numerical solution of D using the spectral method as 
a function of the number of Fourier elements for three different choices of r. 
In the simulation the number of Hermite polynomials Ns was kept constant 
at N, = 30. 



the scope of this paper, we can provide a rough estimation for both methods. 
The WPE-based algorithm involves the solution of two MN x MN system 
of equations (M for the grid size in for the grid size in F). This 

is performed usually in 0((MA^)^)-flopcl, but further examination of the 
matrices involved in the method reveals that the system is sparse and banded, 
reducing the cost of the solutions to 0((MA^)^)-flops. The spectral numerical 
method involves the solution of two sets of recursive systems (one for p and 



one for x) of iV + 1 equations of the form (see Appendix A), 

-(Q„ + Q,-S„+i)"'Qn- 

Although numerically the inverse matrix is never explicitly computed, the 
above operation is numerically equivalent to solving 2M+1 systems of {2M+ 
1) X (2M-I-1) equations, which since all the matrices involved in this equation 



^Using a Gauss-Seidel method, for instance. 
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Figure 4: Error in the numerical solution of D using the WPE method as 
a function of the grid size Ax = l/M^ for three different choices of r. In 
the simulation the number of grid points Np in the F direction was kept 
constant at Np = 30. 



are also sparse and banded, can be done in 0(M^)-flops, so that the total cost 
of the spectral numerical method is 0(A^M^)-flops. The comparison between 
the cost of the numerical methods is done in the following way. By keeping 
the number of A^-elements (either Hermite polynomials or grid points in the 
F-direction) we start with a small number of M-elements (both Fourier and 
X-grid points.) The number of M-elements is then increased such that the 
cost in both numerical methods is increased by (approximately) the same 
factor, i.e., while we double the number of Fourier elements (increasing the 
cost in the spectral algorithm by a factor of eight) we triple the number of grid 
points (increasing the cost in the WPE-based algorithm by a factor of nine). 
In the same manner, while keeping the number of M-elements fixed, while 
we double the number of A^-elements for the WPE-based method (increasing 
the cost by a factor of four) we take 4 N elements for the spectral algorithm 
(increasing the cost also by a factor of four). In Figure [7] and Figure [H] we 
show the results for two different values of r. As we can observe from these 
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Figure 5: Error in the numerical solution of D using the spectral method as 
a function of the number of Hermite polynomials for three different choices 
of r. In the simulation the number of Fourier terms Mg was kept constant 
at Ms = 10. 



results, the convergence of the spectral method is much faster than that of 
the WPE method relative to the numerical cost involved in both methods. 



6. Summary and Discussion. 

We have presented a novel numerical algorithm for computing the effective 
transport properties of the flashing ratchet with continuous Gaussian modu- 
lations (Ornstein-Uhlenbeck process). This numerical algorithm is based on 
a spectral decomposition of the solution to the stationary Fokker-Planck and 
Poisson equations that arise in homogenization theory. The method is shown 
to produce results in agreement with Monte Carlo simulations, with much less 
computational expense. We have also compared this spectral homogenization 
algorithm with a finite volume variation of another computational approach 



due to WPE [17|, Il8|, which can be applied once the continuous modulations 
are discretized into a continuous-time Markov chain. Both algorithms have 
been shown to be theoretically equivalent, and capable of accurately repro- 
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Figure 6: Error in the numerical solution of D using the WPE method as 
a function of the grid size Af ~ 1/Np for three different choices of r. In 
the simulation the number of grid points in the X direction was kept 
constant at = 25. 



ducing the results of Monte Carlo simulations, with the error of our spectral 
method converging to zero more rapidly with increasing computational effort. 
We have also examined to what extent the continuity or discreteness of the 
potential modulations affects the transport properties of the motor particle. 
In one direction, the WPE computational approach is based from the start 
on a discretization of the state space of the random modulations, and we 
have found that with 21 states, the WPE method successfully computes the 
effective drift and diffusivity of the flashing ratchet model over a wide range 
of parameters, except in the white noise limit when the correlation time of 
the modulations is taken small while their amplitude is taken large. Pre- 
sumably a larger number of states are needed for accurate representation by 
the WPE method in this regime. From another perspective, we considered 
how a relatively crude approximation of the continuous Ornstein-Uhlenbeck 
process for the potential modulations in terms of a 2-state (dichotomous) 
Markov chain with the same mean and correlation function affects the trans- 
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Figure 7: Cost comparison between the spectral and the WPE-based (finite- 
volume) numerical algorithms for r = 0.1, ap = 10, and 6 = 0.1 (case 
of enhanced diffusion). Upper panel: The numerical cost of the spectral 
algorithm is estimated as M^, while the numerical cost for the WPE-based 
algorithm is estimated as M^. The number of iV-elements is kept fixed at 
Ns = 20 for the spectral method and Np = 40 for the WPE-based method. 
Bottom panel: The numerical cost of the spectral algorithm is estimated as 
Ng, while the numerical cost for the WPE-based algorithm is estimated as 
Np. The number of M-elements is kept fixed at = 10 for the spectral 
method and ikL = 50 for the WPE-based method. 



port properties of the fiashing ratchet. We found that the dichotomous and 
continuous models produced similar behavior for the motor particle over a 
broad range of parameters, except when either the correlation time or the 
amplitude of the noise is sufficiently large. 

This paper can be extended in several directions, to be explored in future 
work. First, more general types of modulations, not necessarily described by 
Ornstein-Uhlenbeck processes, can be considered. The extension is straight- 
forward for modulations described by reversible diffusions, for which an ap- 
propriate orthonormal basis can be constructed using the eigenfunctions of 
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Figure 8: Cost comparison between the spectral and the WPE-based (finite- 
volume) numerical algorithms for r = 10, = 10, and 6 = 0.1(case of 
suppressed diffusion). Upper panel: The numerical cost of the spectral al- 
gorithm is estimated as M^, while the numerical cost for the WPE-based 
algorithm is estimated as M^. The number of iV-elements is kept fixed at 
Ns = 20 for the spectral method and Np = 40 for the WPE-based method. 
Bottom panel: The numerical cost of the spectral algorithm is estimated as 
Ng, while the numerical cost for the WPE-based algorithm is estimated as 
Np. The number of M-elements is kept fixed at Mg = 10 for the spectral 
method and = 50 for the WPE-based method. 



the generator of the process (e.g. the Hermite polynomials for the Ornstein- 
Uhlenbeck process). Second, we believe that our algorithm can be extended 
to higher dimensional problems and to systems of coupled SDEs/Fokker- 
Planck equations with forcing terms modulated by a stochastic process. This 
would require the use of appropriate tensor products of Hermite polynomi- 
als and Fourier basis functions, together with appropriate preconditioning to 
reduce the computational cost. Third, the rigorous numerical analysis of our 
algorithm, establishing convergence and analyzing its stability properties, is 
a natural next step. 
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Appendix A. The Spectral Numerical Method 

For simplicity in the presentation we choose 0(x) to be = — ^ coswix, 
with Ui = 2iT, although more complex potentials may be considered. The 
spectral representation of (j)'{x) is then simply, 

0'(a;) = 1 (e^^^^ - e'-^-i^) . 

This leads to the following spectral representation of (!23|) in terms of the 
Fourier coefficients {tt^} from Eq. (!25l) . with j = —oo, . . . ,00, 

ap^ujj {irt' - 7it^) - euj]7ri = (A.la) 

I'TZn+i + U7r„ + L+7r„_i = 0, n = 1, 2, . . . (A. lb) 

where 7r„ is an infinite column vector of Fourier coefficients of 7r„ (x), and the 
matrix-vector products above are shorthand for the following operations on 
Fourier coefficients: 



The normalization of the solution, p{x, /) implies furthermore, 

p1 poo p1 poo 

/ / p(x,/)d/dx = / / pFif)J2^nix)H4f)dfdx 
•Jo J —oo Jo J —oo n=0 

oo 

n=0 



< = 1- 
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It can also be noticed from the j = component of equations (lA.lb|) that 
7r° = 0, n = 1, 2, . . .. We now approximate 7r{x) by applying a Galerkin 
truncation to the infinite series at suitable finite values Ns and Mg, 

Ns Ms 
n=0 j=-Ms 

By taking the values of 7r° as given above, the system of equations (123!) 
becomes then a finite system of {2Ms) x (A^^ + 1) linear equations, which can 
be written as, 



Qq TTi + QoTTo = 0, 
Q^7r2 + QlTTi + Qi'^Tro = Bo, 



(A.3a) 



0. 



The matrices Q = {Qi+Ms+ij+Ms+i}, j, I = -Mg, . . . , -1, 1, . . . , then 
take the form, 



[Qo1 



l+Ms + l,j+Ms+l 



-ap-ui ifj = 1 + 1,1 = -Ms, . . . , -2, 1, . . . , M, - 1, 



(TF-coi a J = 1-1,1 = -Ms + 1, . . . , -1, 2, . . . , M, 



otherwise. 



(A.4a) 



[Qo: 



l+Ms+l,j+Ms + l 



-9cof if J = 1,1 = -Ms, ...,-1,1,..., Ms, 
otherwise. 



(A.4b) 



and for n = 1, 2, . . . , Ns, 



[Qn] 



-apVn + 1-ui if J = 1 + 1, I = -Ms, ...,-2,1,..., Ms -1, 



l+Ms+l,j+Ms+l 



apVu + l-oJi if J = 1-1,1 = -Ms + 1, . . . , -1, 2, . . . , M, 







otherwise. 



(A.4c) 
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l+M^+lJ+M^+l 



-opy/n^bJi a j = 1 + 1,1 
apy/n^uji iij = 1-1,1 



-M„ . . . , -2, 1, . . . , M, - 1, 



-M, + 1,...,-1,2,...,M, 







otherwise. 



(A.4d) 



-eojf - nr-^ iij^lj^ -Ms, ...,-l,l,...,M„ 



[Qn. 



l+Ms+l,j+Ms+l 







otherwise. 



(A.4e) 



This system is then solved recursively for ttat^-i, . . . , 7r2 in the form 



with 



TTji — SjiTTn— 1, Tl — Ng, . . . , 1, 



(A.5a) 



(A.5b) 
(A.5c) 



S„ = -(Qn + Q-Sn+i)-^Q+,n = Ng-l,...,2. 
This leaves us with the following set of equations: 

Qo TTi + QoTTo = 0, 

(QrS2 + Ql) TTi + Q+TTo = Bo. 

which can be solved for ttq and tti. Then 7r„, n = 2,3,... can then be 
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recovered using Eq. flA.51) . U is then computed as follows. From (|T7|) . 



p1 poo 

U = - / / 0'(x)/p(a;,/)d/da; 

Jo J-oo 



1 POO 



J-oo 

°° pi POO 



n=0 



pi 

= <P'{x)7Tn{x)aFSi^ndx ^j^q^ 

n=0 

= — (Tf / (j)' {x)Tri{x)dx 
Jo 

fl 1 

since vr"-^ = vr^, where A is the complex conjugate of A G C. 

The solution to the cell problem (fT9|) to compute D is done similarly. We 
express x{Xi f) in its Hermite polynomial decomposition, 



X[x. 



f)=Y,Xn{x)KnU)- 



n=0 



Upon substitution of the above expression in ( |T9l) and using the orthogonality 
of we obtain the following set of equations, 

(t)\x)aFd.j:Xi - Gd^xXo = U, (A.7a) 

4){x)\PiaFdxX2 - QdxxXx + t'^Xx + 0'(a;)(TF9a;Xo = -c^f^', (A.7b) 
- (/:-)kn+i - 4xn - (>C+)lxn-i = 0, n = 2, 3, . . . (A.7c) 
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where the operators £ and are defined in Eq. (^^, and f denotes the 
adjoint. By decomposing Xn{x) in terms of its Fourier series, 



Xn[X) 



the set of equations f lA.7p becomes the hnear system, for j = — oo. 



:n\x2 + Lixi + (L+)ixo 



-U(5j,o, 

0, n = 2,3, ... 



where x„ denotes the sequence of Fourier coefficients of Xn{x), and L~,L, and 
L"*" are defined in Eq. (1A.2I) . In order to solve this infinite set of equations, 
we apply a Galerkin truncation to the spectral decomposition of /) at 
some appropriate values Ms and Ng: 

Ns Ms 
n=0 j=-Ms 

SO we get the finite {2Ms + l)(A^s + 1) set of algebraic equations. 



(Q;)^Xn+l + QlXn + (Qt^Xn-l 



-Bo, 

0, n = 2,3. 



AT, 



(A.8a) 
(A.8b) 
(A.8c) 



where the matrices Q„, Q„, and were defined in flA.4p . and 



0\l+Ah + l 



-U5;.o, 



l+Ms+l 



^A 



The supermatrix implicitly defined by the left hand side of this system does 
have a zero eigenvalue, with right eigenvector {do, 0, . . . , 0)''' and left eigen- 
vector (ttq, TTi, . . . , tttvJ, all inherited from discretization of the operator C 
Solvability of f lA.SP follows from verifying that the supervector composed of 
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the right hand sides is indeed orthogonal to the null left eigenvector of the 
supermatrix, i.e., 

TTo ■ So + TTi ■ Si = 

as follows from the definitions of these vectors and the formula flA.6|) for the 
effective velocity U. A unique solution is obtained by imposing Xo = 0, the 
analogue of the constraint f l20|) for the continuous formulation. 

The system of equations is then solved similarly as we did for 7r„(x), 
beginning with writing: 



'X.n 



n 



Ns, - 1, 



,2, 



.n = N, 



(A.9a) 
(A.9b) 

(A.9c) 
(A.IO) 



From the equation for = 1, we can write, 

Xi = Bi + ZiXo, 

where Zi is defined by the same formula as for n > 2 in Eq. (lA.9cp . and 

Bi= (PrZ2 + Pi)"'Bi, 
Substituting this expression in the equation for n = yields 

ZoXo = ^0, (A.ll) 

with 

Zo = Pq Zi + Po, 

and 

Bo = Bo — Po Bi. 

We handle the degeneracy of the matrix Zo by imposing Xo = ^"^^ solving 
for the remaining components of Xo- '^^^ solution is completed by computing 
{Xn}n=i through Eqs. f OTTOj) and HKMl . 

From Eq. f l?T]) . a similar calculation to that in Eq. flA.6p yields 



Ms 



D = e + i—^V^^ 

n=0 

Ns Ms 

+47rei^ jxIttL 

71=0 j = -Ms 



Ej j+1 j J— 1 I j j+1 j j—1 



\J=-Ms 



(A.12) 
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Appendix B. Discrete-State Approximation of the Ornstein-Uhlenbeck 
Process 



In our numerical experiments, we discretize the state space for the noise 
variable with Np = 2np + l states, equally spaced with interval A/, and cen- 
tered about 0: = {-n^A/, -{np - 1) A/, . . . , 0, . . . , (rij. - 1) A/, npAf}. 
Equivalently, 5*1;. = {fn}n=i with fn = {n — np — 1)A/. We will approximate 
f l57|) at the grid points /„ by first approximating the derivative by a centered 
finite difference at the points / = /« + A//2 and / = /« — A//2, 

&t ^ T Af 
where. 

Unit) = Uifn,t) 

is just the point evaluation of the function uif,t) at the grid point / = /„. 
Next, we approximate the derivative at / = /„ ± A//2 once again by a 
centered difference this time around the grid points / = /„± A/ and / = /„. 
The final approximation can be written as, 

~ Kn,n+lUn+lit) — (-ft'n,n+l + Kn,n-l) Unit) + Kn,n-lUn-lit) , 

with, 

LC , - _^p-/3(V(/n+A//2)-V(/„)) (-n 1 N 



K , _ „-^(V(/n-A//2)-y(/„)) 



The approximation of the backward-Kolmogorov equation can be expressed 
as, 

^n(t) = Luit), 

where the entries of the matrix L are given as 

Kn,n+i if n' = n + 1,1 < n,n' < Np — 1 

Kn,n-i a n' = n — 1,2 < n,n' < Np 

-Kn,n+l - Kn,n-1 if u' = U, 2 < U < Np - 1 

—Ki 2 a n' = n = 1 

-Knp,Np-i iin' = n = Np 

otherwise. 



[L]n,n' 
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which, for any choice of A/ > 0, defines a Markov jump process with space 
state defined by the grid points {fn}n=i ^^cl jump rates between /„ and /„±i 
given by Kn,n±i- Moreover, it can be easily checked that the vector, 

N„ = e-'^^(^"), 

solves the equation, 

L*7r = 0, 

where L* is the adjoint matrix of L. This means that tt is an invariant dis- 
tribution of the Markov jump process and is consistent with the invariant 
distribution of the continuous process F{t). Moreover, it can also be shown 
that the Markov jump process satisfies the detailed balance condition with 
respect to this invariant measure, also consistently with the continuous pro- 
cess. 

We must next truncate the infinite state space, and impose boundary condi- 
tions. Since the OU-process has a stationary Gaussian distribution p{f) ~ 
e'^'f we choose the last grid point to be such that, 

where 5 is a small number. In practice we choose S = 10"^'^. The discrete- 
state, Markov jump approximation of the process F^(t) has state space Sp = 
{—npAf, —{rip — 1)A/, . . . , 0, . . . , {rip — 1)A/, np/S.f} and transition rates 
given by fIB.ip . 
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